# data vizlibrary(ggplot2)library(ggh4x)library(ggOceanMaps)library(patchwork)library(viridis)library(ggpubr)library(grid)library(viridis)library(ggnewscale)# tablelibrary(gt)# maplibrary(rnaturalearth)library(ggmap)library(ggsn)library(sf)library(sp)library(smoothr)# kernel densitylibrary(eks)library(ggdensity)# statlibrary(Hmisc)library(circular)library(CircStats)# charlibrary(stringr)# solar anglelibrary(oce)# data wranglinglibrary(tidyr)library(dplyr)library(purrr)library(forcats)library(lubridate)# modellibrary(lme4)library(fitdistrplus)
2 Setting up custom function
2.1windrose
Code
windrose <-function(data_to_plot,grid =NULL,set_title =NULL,legend_position ="none",facet = F) {# this code comes from Roxanne uniqhours <-1:24* (360/24)# trick to align hours om the graph data_to_plot <-rbind(data_to_plot[2:nrow(data_to_plot), ], data_to_plot[1, ])for (i in1:24) {# turn hours to radiansif (i ==1) { temp <-rep(uniqhours[i], data_to_plot$nb_ind_hour[i]) day_night <-rep("night", data_to_plot$nb_ind_hour[i]) } else { temp <-c(temp, rep(uniqhours[i], data_to_plot$nb_ind_hour[i])) day_night <-c(day_night, rep(if_else(between(i, 7, 20), "day", "night"), data_to_plot$nb_ind_hour[i] )) } } data2 <-data.frame(direction = temp) deg <-15# choose bin size (degrees/bin) dir.breaks <-seq(0- (deg /2), 360+ (deg /2), deg) # define the range of each bin# assign each direction to a bin range dir.binned <-cut(data2$direction,breaks = dir.breaks,ordered_result =TRUE )# generate pretty lables dir.labels <-as.character(c(seq(0, 360- deg, by = deg), 0))# replace ranges with pretty bin lableslevels(dir.binned) <- dir.labels# Assign bin names to the original data set data2$dir.binned <- dir.binned# add origin if anyif (facet) { data2$origin <-unique(data_to_plot$origin) }# set up max value maxvalue <-35# initialise the plot plt.dirrose_2 <-ggplot()# check if gridif (!is.null(grid)) { plt.dirrose_2 <- plt.dirrose_2 +geom_hline(yintercept = grid,colour ="grey20",linewidth = .2 ) } plt.dirrose_2 <- plt.dirrose_2 +geom_vline(xintercept =c(seq(1, 24, 2)),colour ="grey30",linewidth =0.2 ) +# 24 vertical lines at center of the 30? ranges.geom_hline(yintercept = maxvalue,colour ="white",linewidth = .5 ) +# Darker horizontal line as the top border (max).# On top of everything we place the histogram bars.geom_bar(data = data2,aes(x = dir.binned, fill = day_night),width =1,colour ="black",linewidth =0.3 ) +# geom_bar(data = data2, aes(x = dir.binned2), width = 1, colour="black", size = 0.3,fill="salmon",alpha=0.9) +scale_x_discrete(drop =FALSE,labels =c(0, "", 2, "", 4, "", 6, "", 8, "", 10, "", 12, "", 14, "", 16, "", 18, "", 20, "", 22, "" ) ) +scale_fill_manual(values =c("white", "darkgrey"), name ="Time of day") +labs(x ="Time (hours)", y ="Count", title = set_title) +coord_polar(start =-(deg /2) * (pi /180))# if facetif (facet) { plt.dirrose_2 <- plt.dirrose_2 +facet_wrap2(. ~ origin) }# Wraps the histogram into a windrose plt.dirrose_2 <- plt.dirrose_2 +theme_bw() +theme(legend.position = legend_position,panel.grid.minor =element_blank(),panel.grid.major =element_blank(),panel.background =element_blank(),legend.key.size =unit(0.5, "cm") )# returnreturn(plt.dirrose_2) }
getNightLength <-function(survey_date, latitude, longitude) {# survey_date <- strptime(survey_date, "%d/%m/%y") time_origin <-strptime("01/01/2000 12:00", "%d/%m/%Y %H:%M") n <-as.numeric(difftime(survey_date, time_origin, units ="days")) Jstar <- n - longitude /360 M <- (357.5291+0.98560028* Jstar) %%360 C <-1.9148*sin(M *2* pi /360) +0.0200*sin(2* M *2* pi /360) +0.0003*sin(3* M *2* pi /360) lambda <- (M + C +180+102.9372) %%360 Jtransit <-as.double(2451545.000) +as.double(Jstar) +as.double(0.0053*sin(M *2* pi /360)) -as.double(0.0069*sin(lambda *4* pi /360)) sindelta <-sin(lambda *2* pi /360) *sin(23.44*2* pi /360) delta <-asin(sindelta) *360/ (2* pi) cos_omega <- (sin(-0.83*2* pi /360) - (sindelta *sin(latitude *2* pi /360))) / (cos(latitude *2* pi /360) *cos(delta *2* pi /360)) omega <-acos(cos_omega) *360/ (2* pi) Jrise <- Jtransit - omega /360 Jset <- Jtransit + omega /360return((1- (Jset - Jrise)) *24)}
3 Import data
Let’s load data_dive, i.e. the output of data_wrangling.qmd, and filter only on animals leaving from Ano Nuevo.
Code
# import the datadata_dive <-readRDS("../export/data_dive.rds")# filter on seals departing from ANNUdata_dive <- data_dive %>%filter(DepartureLocation =="ANNU")
Code
# get solar angle when data locationsolar_angle_inter <- data_dive %>%filter(!is.na(Lat)) %>%mutate(solar_angle =sunAngle(date, Long, Lat)$altitude) %>%select(id, DiveNumber, solar_angle)# merge that with data_divedata_dive <-merge(data_dive, solar_angle_inter,by =c("id", "DiveNumber"),all.x = T)# add day-night# https://sciencepickle.com/earth-systems/sun-earth-connection/declination-circles/sunrise-sunset-and-twilight/data_dive <- data_dive %>%mutate(day_night =if_else(solar_angle <-18, "night", "day"))
Let’s add a condition based on the difference between the maximum depth reached and the bathymetry to define a benthic dive, i.e. if this difference is above xxx m then is not considered as a benthic dive.
Code
data_dive <- data_dive %>%# change DiveType to include condition on difference# between max depth and bathymetry (<50m)mutate(DiveType_50 = data.table::fifelse(DiveType ==3& (abs(bathy) - Maxdepth) >50, 0, DiveType),DiveTypeName_50 = data.table::fifelse( DiveTypeName =="Benthic"& (abs(bathy) - Maxdepth) >50,"Transit", DiveTypeName ),# same with <100mDiveType_100 = data.table::fifelse(DiveType ==3& (abs(bathy) - Maxdepth) >100, 0, DiveType),DiveTypeName_100 = data.table::fifelse( DiveTypeName =="Benthic"& (abs(bathy) - Maxdepth) >100,"Transit", DiveTypeName ),# same with <150mDiveType_150 = data.table::fifelse(DiveType ==3& (abs(bathy) - Maxdepth) >150, 0, DiveType),DiveTypeName_150 = data.table::fifelse( DiveTypeName =="Benthic"& (abs(bathy) - Maxdepth) >150,"Transit", DiveTypeName ),# same with <200mDiveType_200 = data.table::fifelse(DiveType ==3& (abs(bathy) - Maxdepth) >200, 0, DiveType),DiveTypeName_200 = data.table::fifelse( DiveTypeName =="Benthic"& (abs(bathy) - Maxdepth) >200,"Transit", DiveTypeName ) )
4 Maps
Code
# check if background_ggoceanmaps existif (file.exists("../export/background_ggoceanmap.rds")) { trip <-readRDS("../export/background_ggoceanmap.rds")} else {# using ggOceanMaps trip <-basemap(limits =c(170, -110, 20, 59),bathymetry =TRUE,shapefiles ="Arctic",rotate =TRUE,grid.col =NA )# Make the graticules: lims <-attributes(trip)$limits graticule <- sf::st_graticule(c(lims[1], lims[3], lims[2], lims[4]),crs =attributes(trip)$proj,lon =seq(-180, 180, 45),lat =seq(-90, 90, 10) )# Plot trip <- trip +geom_sf(data = graticule, color ="grey50") trip$layers <- trip$layers[c(1, 3, 2)]# save resultsaveRDS(trip, "../export/background_ggoceanmap.rds")}
# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat) & DiveTypeName =="Benthic") %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName)
Figure 1: Kernel density plots of the all benthic dives performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.
Code
# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat) & DiveTypeName_50 =="Benthic") %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName)
Code
# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plottest <- trip +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# titlelabs(title ="Benthic Dives with bathy-max_depth<50",subtitle =paste(nrow(df_kernel_dens), "dives") )# reorder layerstest$layers <- test$layers[c(1, 2, 4, 3)]# plottest
Figure 2: Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 50 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.
Code
# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat) & DiveTypeName_100 =="Benthic") %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName)
Code
# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plottest <- trip +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# titlelabs(title ="Benthic Dives with bathy-max_depth<100",subtitle =paste(nrow(df_kernel_dens), "dives") )# reorder layerstest$layers <- test$layers[c(1, 2, 4, 3)]# plottest
Figure 3: Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 100 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.
Code
# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat) & DiveTypeName_150 =="Benthic") %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName)
Code
# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plottest <- trip +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# titlelabs(title ="Benthic Dives with bathy-max_depth<150",subtitle =paste(nrow(df_kernel_dens), "dives") )# reorder layerstest$layers <- test$layers[c(1, 2, 4, 3)]# plottest
Figure 4: Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 150 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.
Code
# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat) & DiveTypeName_200 =="Benthic") %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName)
Code
# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plottest <- trip +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# titlelabs(title ="Benthic Dives with bathy-max_depth<200",subtitle =paste(nrow(df_kernel_dens), "dives") )# reorder layerstest$layers <- test$layers[c(1, 2, 4, 3)]# plottest
Figure 5: Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 200 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.
# let's add a column with the local_timedata_windrose <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName =="Benthic")# figurewindrose_benthic <- data_windrose %>%group_by(id) %>%filter(DiveTypeName =="Benthic") %>%mutate(time =as.numeric(format(date_tz, format ="%H"))) %>%group_by(time) %>%summarise(nb_ind_hour =n()) %>%mutate(origin ="Benthic") %>%windrose(., grid =c(3000, 6000, 9000), legend_position ="top", facet = T)# displaywindrose_benthic
Figure 6: Circular histogram plots displaying the time when they perform benthic dives across their foraging trip.
Code
# perform testtest_benthic <-r.test(conversion.circular(circular( data_windrose %>%group_by(id) %>%filter(DiveTypeName =="Benthic") %>%mutate(time =hour(date_tz) +minute(date_tz) /60+second(date_tz) / (60*60)) %>%pull(time),units ="hours" ),units ="radians"))# convert back into human reading hourconv_test <-conversion.circular(circular(test_benthic$r.bar,units ="rad" ),units ="hours")# printprint(paste0("The average time is ", round(conv_test[[1]], 1), " (p-value: ", round(test_benthic$p.value, 1), ")"))
[1] "The average time is 0.5 (p-value: 0)"
Code
# let's add a column with the local_timedata_windrose <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_50 =="Benthic")# figurewindrose_benthic <- data_windrose %>%group_by(id) %>%filter(DiveTypeName =="Benthic") %>%mutate(time =as.numeric(format(date_tz, format ="%H"))) %>%group_by(time) %>%summarise(nb_ind_hour =n()) %>%mutate(origin ="Benthic") %>%windrose(., grid =c(1000, 2000, 3000), legend_position ="top", facet = T)# displaywindrose_benthic
Figure 7: Circular histogram plots displaying the time when they perform benthic dives (diff<50) across their foraging trip.
Code
# perform testtest_benthic <-r.test(conversion.circular(circular( data_windrose %>%group_by(id) %>%filter(DiveTypeName_50 =="Benthic") %>%mutate(time =hour(date_tz) +minute(date_tz) /60+second(date_tz) / (60*60)) %>%pull(time),units ="hours" ),units ="radians"))# convert back into human reading hourconv_test <-conversion.circular(circular(test_benthic$r.bar,units ="rad" ),units ="hours")# printprint(paste0("The average time is ", round(conv_test[[1]], 1), " (p-value: ", round(test_benthic$p.value, 1), ")"))
[1] "The average time is 0.1 (p-value: 0)"
Code
# let's add a column with the local_timedata_windrose <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_100 =="Benthic")# figurewindrose_benthic <- data_windrose %>%group_by(id) %>%filter(DiveTypeName =="Benthic") %>%mutate(time =as.numeric(format(date_tz, format ="%H"))) %>%group_by(time) %>%summarise(nb_ind_hour =n()) %>%mutate(origin ="Benthic") %>%windrose(., grid =c(1000, 2000, 3000), legend_position ="top", facet = T)# displaywindrose_benthic
Figure 8: Circular histogram plots displaying the time when they perform benthic dives (diff<100) across their foraging trip.
Code
# perform testtest_benthic <-r.test(conversion.circular(circular( data_windrose %>%group_by(id) %>%filter(DiveTypeName_100 =="Benthic") %>%mutate(time =hour(date_tz) +minute(date_tz) /60+second(date_tz) / (60*60)) %>%pull(time),units ="hours" ),units ="radians"))# convert back into human reading hourconv_test <-conversion.circular(circular(test_benthic$r.bar,units ="rad" ),units ="hours")# printprint(paste0("The average time is ", round(conv_test[[1]], 1), " (p-value: ", round(test_benthic$p.value, 1), ")"))
[1] "The average time is 0.1 (p-value: 0)"
Code
# let's add a column with the local_timedata_windrose <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_150 =="Benthic")# figurewindrose_benthic <- data_windrose %>%group_by(id) %>%filter(DiveTypeName =="Benthic") %>%mutate(time =as.numeric(format(date_tz, format ="%H"))) %>%group_by(time) %>%summarise(nb_ind_hour =n()) %>%mutate(origin ="Benthic") %>%windrose(., grid =c(1000, 2000, 3000), legend_position ="top", facet = T)# displaywindrose_benthic
Figure 9: Circular histogram plots displaying the time when they perform benthic dives (diff<150) across their foraging trip.
Code
# perform testtest_benthic <-r.test(conversion.circular(circular( data_windrose %>%group_by(id) %>%filter(DiveTypeName_150 =="Benthic") %>%mutate(time =hour(date_tz) +minute(date_tz) /60+second(date_tz) / (60*60)) %>%pull(time),units ="hours" ),units ="radians"))# convert back into human reading hourconv_test <-conversion.circular(circular(test_benthic$r.bar,units ="rad" ),units ="hours")# printprint(paste0("The average time is ", round(conv_test[[1]], 1), " (p-value: ", round(test_benthic$p.value, 1), ")"))
[1] "The average time is 0.1 (p-value: 0)"
Code
# let's add a column with the local_timedata_windrose <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_200 =="Benthic")# figurewindrose_benthic <- data_windrose %>%group_by(id) %>%filter(DiveTypeName =="Benthic") %>%mutate(time =as.numeric(format(date_tz, format ="%H"))) %>%group_by(time) %>%summarise(nb_ind_hour =n()) %>%mutate(origin ="Benthic") %>%windrose(., grid =c(1000, 2000, 3000), legend_position ="top", facet = T)# displaywindrose_benthic
Figure 10: Circular histogram plots displaying the time when they perform benthic dives (diff<200) across their foraging trip.
Code
# perform testtest_benthic <-r.test(conversion.circular(circular( data_windrose %>%group_by(id) %>%filter(DiveTypeName_200 =="Benthic") %>%mutate(time =hour(date_tz) +minute(date_tz) /60+second(date_tz) / (60*60)) %>%pull(time),units ="hours" ),units ="radians"))# convert back into human reading hourconv_test <-conversion.circular(circular(test_benthic$r.bar,units ="rad" ),units ="hours")# printprint(paste0("The average time is ", round(conv_test[[1]], 1), " (p-value: ", round(test_benthic$p.value, 1), ")"))
[1] "The average time is 0.1 (p-value: 0)"
6 Boxplot of the benthic dive rate (day vs. night)
Since we thought of accoutering for the variation in dive duration to explore how benthic dive rate varies across day-time vs night-time, I just looked at the difference in dive duration between both period of the day.
Code
data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%ggplot(aes(x = day_night, y = Dduration)) +# the boxplotgeom_boxplot(outlier.shape =NA) +# zoomcoord_cartesian(ylim =c(0, 2600)) +# labelslabs(x ="Time of the day", y ="Dive duration (s)")
Code
broom::tidy(wilcox.test( data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# at nightfilter(day_night =="night") %>%# get dive durationpull(Dduration), data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# at nightfilter(day_night =="day") %>%# get dive durationpull(Dduration)))
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName =="Benthic") %>%mutate(hour = lubridate::hour(date),day =as.Date(date) ) %>%group_by(id, day, hour, day_night) %>%summarise(nb_benthic_dive =n(),mean_dive_duration =mean(Dduration),.groups ="drop" ) %>%group_by(id, day_night) %>%summarise(nb_average_benthic_dive =mean(nb_benthic_dive),mean_dive_duration =mean(mean_dive_duration),.groups ="drop" )data_plot %>%ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +geom_boxplot(outlier.shape =NA) +labs(x ="Time of the day",y ="Average number of benthic dives per seal per hour" )
Figure 11: Boxplot of the number of (all) benthic dives per hour and individual accross day-time and night-time
Model
Code
# determine what distribution would best describe our datadescdist(data_plot$nb_average_benthic_dive, discrete =FALSE)
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName =="Benthic") %>%# I'm not so sure about this part..mutate(period_duration =if_else( day_night =="night",getNightLength(date, Lat, Long),24-getNightLength(date, Lat, Long) )) %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(),mean_dive_duration =mean(Dduration),mean_period_duration =mean(period_duration),.groups ="drop" )
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_50 =="Benthic") %>%mutate(hour = lubridate::hour(date),day =as.Date(date) ) %>%group_by(id, day, hour, day_night) %>%summarise(nb_benthic_dive =n(),mean_dive_duration =mean(Dduration),.groups ="drop" ) %>%group_by(id, day_night) %>%summarise(nb_average_benthic_dive =mean(nb_benthic_dive),mean_dive_duration =mean(mean_dive_duration),.groups ="drop" )data_plot %>%ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +geom_boxplot(outlier.shape =NA) +labs(x ="Time of the day",y ="Average number of benthic dives per seal per hour" )
Figure 13: Boxplot of the number of (diff<50) benthic dives per hour and individual accross day-time and night-time
Model
Code
# determine what distribution would best describe our datadescdist(data_plot$nb_average_benthic_dive, discrete =FALSE)
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_50 =="Benthic") %>%# I'm not so sure about this part..mutate(period_duration =if_else( day_night =="night",getNightLength(date, Lat, Long),24-getNightLength(date, Lat, Long) )) %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(),mean_dive_duration =mean(Dduration),mean_period_duration =mean(period_duration),.groups ="drop" )
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_100 =="Benthic") %>%mutate(hour = lubridate::hour(date),day =as.Date(date) ) %>%group_by(id, day, hour, day_night) %>%summarise(nb_benthic_dive =n(),mean_dive_duration =mean(Dduration),.groups ="drop" ) %>%group_by(id, day_night) %>%summarise(nb_average_benthic_dive =mean(nb_benthic_dive),mean_dive_duration =mean(mean_dive_duration),.groups ="drop" )data_plot %>%ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +geom_boxplot(outlier.shape =NA) +labs(x ="Time of the day",y ="Average number of benthic dives per seal per hour" )
Figure 15: Boxplot of the number of (diff<100) benthic dives per hour and individual accross day-time and night-time
Model
Code
# determine what distribution would best describe our datadescdist(data_plot$nb_average_benthic_dive, discrete =FALSE)
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_100 =="Benthic") %>%# I'm not so sure about this part..mutate(period_duration =if_else( day_night =="night",getNightLength(date, Lat, Long),24-getNightLength(date, Lat, Long) )) %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(),mean_dive_duration =mean(Dduration),mean_period_duration =mean(period_duration),.groups ="drop" )
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_150 =="Benthic") %>%mutate(hour = lubridate::hour(date),day =as.Date(date) ) %>%group_by(id, day, hour, day_night) %>%summarise(nb_benthic_dive =n(),mean_dive_duration =mean(Dduration),.groups ="drop" ) %>%group_by(id, day_night) %>%summarise(nb_average_benthic_dive =mean(nb_benthic_dive),mean_dive_duration =mean(mean_dive_duration),.groups ="drop" )data_plot %>%ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +geom_boxplot(outlier.shape =NA) +labs(x ="Time of the day",y ="Average number of benthic dives per seal per hour" )
Figure 17: Boxplot of the number of (diff<150) benthic dives per hour and individual accross day-time and night-time
Model
Code
# determine what distribution would best describe our datadescdist(data_plot$nb_average_benthic_dive, discrete =FALSE)
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_150 =="Benthic") %>%# I'm not so sure about this part..mutate(period_duration =if_else( day_night =="night",getNightLength(date, Lat, Long),24-getNightLength(date, Lat, Long) )) %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(),mean_dive_duration =mean(Dduration),mean_period_duration =mean(period_duration),.groups ="drop" )
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_200 =="Benthic") %>%mutate(hour = lubridate::hour(date),day =as.Date(date) ) %>%group_by(id, day, hour, day_night) %>%summarise(nb_benthic_dive =n(),mean_dive_duration =mean(Dduration),.groups ="drop" ) %>%group_by(id, day_night) %>%summarise(nb_average_benthic_dive =mean(nb_benthic_dive),mean_dive_duration =mean(mean_dive_duration),.groups ="drop" )data_plot %>%ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +geom_boxplot(outlier.shape =NA) +labs(x ="Time of the day",y ="Average number of benthic dives per seal per hour" )
Figure 19: Boxplot of the number of (diff<200) benthic dives per hour and individual accross day-time and night-time
Model
Code
# determine what distribution would best describe our datadescdist(data_plot$nb_average_benthic_dive, discrete =FALSE)
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_200 =="Benthic") %>%# I'm not so sure about this part..mutate(period_duration =if_else( day_night =="night",getNightLength(date, Lat, Long),24-getNightLength(date, Lat, Long) )) %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(),mean_dive_duration =mean(Dduration),mean_period_duration =mean(period_duration),.groups ="drop" )
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName =="Benthic") %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_benthic_dive)) +geom_boxplot(outlier.shape =NA) +coord_cartesian(ylim =c(0, 1000)) +labs(x ="Time of the day",y ="Number of benthic dives per individual" )
Figure 21: Boxplot of the average number of (all) benthic dives per hour and individual accross day-time and night-time
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_50 =="Benthic") %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_benthic_dive)) +geom_boxplot(outlier.shape =NA) +coord_cartesian(ylim =c(0, 250)) +labs(x ="Time of the day",y ="Number of benthic dives per individual" )
Figure 22: Boxplot of the average number of (diff<50) benthic dives per hour and individual accross day-time and night-time
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_100 =="Benthic") %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_benthic_dive)) +geom_boxplot(outlier.shape =NA) +coord_cartesian(ylim =c(0, 250)) +labs(x ="Time of the day",y ="Number of benthic dives per individual" )
Figure 23: Boxplot of the average number of (diff<100) benthic dives per hour and individual accross day-time and night-time
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_150 =="Benthic") %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_benthic_dive)) +geom_boxplot(outlier.shape =NA) +coord_cartesian(ylim =c(0, 250)) +labs(x ="Time of the day",y ="Number of benthic dives per individual" )
Figure 24: Boxplot of the average number of (diff<150) benthic dives per hour and individual accross day-time and night-time
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_200 =="Benthic") %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_benthic_dive)) +geom_boxplot(outlier.shape =NA) +coord_cartesian(ylim =c(0, 250)) +labs(x ="Time of the day",y ="Number of benthic dives per individual" )
Figure 25: Boxplot of the average number of (diff<200) benthic dives per hour and individual accross day-time and night-time
---title: "Benthic Dives Investigations"author: - name: "Astarte Brown" - name: "Conner Hale" - name: "Joffrey Joumaa"date: "`r invisible(Sys.setlocale(locale = 'C')); format(Sys.Date(), format = '%B %d, %Y')`"format: html: toc: true toc-location: left number-sections: true smooth-scroll: true code-fold: true code-tools: true df-print: paged fig-align: "center" highlight-style: arrow self-contained: trueexecute: echo: true cache: false warning: falsetheme: light: flatly dark: darklyknitr: opts_chunk: message: false rownames.print: false tidy: styler---# Import library```{r}# data vizlibrary(ggplot2)library(ggh4x)library(ggOceanMaps)library(patchwork)library(viridis)library(ggpubr)library(grid)library(viridis)library(ggnewscale)# tablelibrary(gt)# maplibrary(rnaturalearth)library(ggmap)library(ggsn)library(sf)library(sp)library(smoothr)# kernel densitylibrary(eks)library(ggdensity)# statlibrary(Hmisc)library(circular)library(CircStats)# charlibrary(stringr)# solar anglelibrary(oce)# data wranglinglibrary(tidyr)library(dplyr)library(purrr)library(forcats)library(lubridate)# modellibrary(lme4)library(fitdistrplus)```# Setting up custom function## `windrose````{r}windrose <-function(data_to_plot,grid =NULL,set_title =NULL,legend_position ="none",facet = F) {# this code comes from Roxanne uniqhours <-1:24* (360/24)# trick to align hours om the graph data_to_plot <-rbind(data_to_plot[2:nrow(data_to_plot), ], data_to_plot[1, ])for (i in1:24) {# turn hours to radiansif (i ==1) { temp <-rep(uniqhours[i], data_to_plot$nb_ind_hour[i]) day_night <-rep("night", data_to_plot$nb_ind_hour[i]) } else { temp <-c(temp, rep(uniqhours[i], data_to_plot$nb_ind_hour[i])) day_night <-c(day_night, rep(if_else(between(i, 7, 20), "day", "night"), data_to_plot$nb_ind_hour[i] )) } } data2 <-data.frame(direction = temp) deg <-15# choose bin size (degrees/bin) dir.breaks <-seq(0- (deg /2), 360+ (deg /2), deg) # define the range of each bin# assign each direction to a bin range dir.binned <-cut(data2$direction,breaks = dir.breaks,ordered_result =TRUE )# generate pretty lables dir.labels <-as.character(c(seq(0, 360- deg, by = deg), 0))# replace ranges with pretty bin lableslevels(dir.binned) <- dir.labels# Assign bin names to the original data set data2$dir.binned <- dir.binned# add origin if anyif (facet) { data2$origin <-unique(data_to_plot$origin) }# set up max value maxvalue <-35# initialise the plot plt.dirrose_2 <-ggplot()# check if gridif (!is.null(grid)) { plt.dirrose_2 <- plt.dirrose_2 +geom_hline(yintercept = grid,colour ="grey20",linewidth = .2 ) } plt.dirrose_2 <- plt.dirrose_2 +geom_vline(xintercept =c(seq(1, 24, 2)),colour ="grey30",linewidth =0.2 ) +# 24 vertical lines at center of the 30? ranges.geom_hline(yintercept = maxvalue,colour ="white",linewidth = .5 ) +# Darker horizontal line as the top border (max).# On top of everything we place the histogram bars.geom_bar(data = data2,aes(x = dir.binned, fill = day_night),width =1,colour ="black",linewidth =0.3 ) +# geom_bar(data = data2, aes(x = dir.binned2), width = 1, colour="black", size = 0.3,fill="salmon",alpha=0.9) +scale_x_discrete(drop =FALSE,labels =c(0, "", 2, "", 4, "", 6, "", 8, "", 10, "", 12, "", 14, "", 16, "", 18, "", 20, "", 22, "" ) ) +scale_fill_manual(values =c("white", "darkgrey"), name ="Time of day") +labs(x ="Time (hours)", y ="Count", title = set_title) +coord_polar(start =-(deg /2) * (pi /180))# if facetif (facet) { plt.dirrose_2 <- plt.dirrose_2 +facet_wrap2(. ~ origin) }# Wraps the histogram into a windrose plt.dirrose_2 <- plt.dirrose_2 +theme_bw() +theme(legend.position = legend_position,panel.grid.minor =element_blank(),panel.grid.major =element_blank(),panel.background =element_blank(),legend.key.size =unit(0.5, "cm") )# returnreturn(plt.dirrose_2) }```## `getNightLength`Based on here: <https://stackoverflow.com/questions/59566553/calculate-night-length-in-r>```{r}getNightLength <-function(survey_date, latitude, longitude){# survey_date <- strptime(survey_date, "%d/%m/%y") time_origin <-strptime("01/01/2000 12:00", "%d/%m/%Y %H:%M") n <-as.numeric(difftime(survey_date, time_origin, units ="days")) Jstar <- n - longitude /360 M <- (357.5291+0.98560028* Jstar) %%360 C <-1.9148*sin(M *2*pi/360) +0.0200*sin(2*M *2*pi/360) +0.0003*sin(3*M *2*pi/360) lambda <- (M + C +180+102.9372) %%360 Jtransit <-as.double(2451545.000) +as.double(Jstar) +as.double(0.0053*sin(M *2*pi/360)) -as.double(0.0069*sin(lambda *4*pi/360)) sindelta <-sin(lambda *2*pi/360) *sin(23.44*2*pi/360) delta <-asin(sindelta) *360/(2*pi) cos_omega <- (sin(-0.83*2*pi/360) - (sindelta *sin(latitude *2*pi/360)))/ (cos(latitude *2*pi/360) *cos(delta *2*pi/360)) omega <-acos(cos_omega) *360/ (2*pi) Jrise <- Jtransit - omega /360 Jset <- Jtransit + omega /360return((1- (Jset - Jrise)) *24)}```# Import dataLet's load `data_dive`, *i.e.* the output of `data_wrangling.qmd`, and filter only on animals leaving from Ano Nuevo.```{r}# import the datadata_dive <-readRDS("../export/data_dive.rds")# filter on seals departing from ANNUdata_dive <- data_dive %>%filter(DepartureLocation =="ANNU")``````{r}# get solar angle when data locationsolar_angle_inter <- data_dive %>%filter(!is.na(Lat)) %>%mutate(solar_angle =sunAngle(date,Long,Lat)$altitude) %>%select(id, DiveNumber, solar_angle)# merge that with data_divedata_dive =merge(data_dive, solar_angle_inter,by =c("id", "DiveNumber"), all.x = T)# add day-night# https://sciencepickle.com/earth-systems/sun-earth-connection/declination-circles/sunrise-sunset-and-twilight/data_dive <- data_dive %>%mutate(day_night =if_else(solar_angle <-18, "night", "day"))```Let's add a condition based on the difference between the maximum depth reached and the bathymetry to define a benthic dive, *i.e.* if this difference is above xxx m then is not considered as a benthic dive.```{r}data_dive = data_dive %>%# change DiveType to include condition on difference# between max depth and bathymetry (<50m)mutate(DiveType_50 = data.table::fifelse(DiveType ==3& (abs(bathy) - Maxdepth) >50, 0, DiveType),DiveTypeName_50 = data.table::fifelse( DiveTypeName =="Benthic"& (abs(bathy) - Maxdepth) >50,"Transit", DiveTypeName ),# same with <100mDiveType_100 = data.table::fifelse(DiveType ==3& (abs(bathy) - Maxdepth) >100, 0, DiveType),DiveTypeName_100 = data.table::fifelse( DiveTypeName =="Benthic"& (abs(bathy) - Maxdepth) >100,"Transit", DiveTypeName ),# same with <150mDiveType_150 = data.table::fifelse(DiveType ==3& (abs(bathy) - Maxdepth) >150, 0, DiveType),DiveTypeName_150 = data.table::fifelse( DiveTypeName =="Benthic"& (abs(bathy) - Maxdepth) >150,"Transit", DiveTypeName ),# same with <200mDiveType_200 = data.table::fifelse(DiveType ==3& (abs(bathy) - Maxdepth) >200, 0, DiveType),DiveTypeName_200 = data.table::fifelse( DiveTypeName =="Benthic"& (abs(bathy) - Maxdepth) >200,"Transit", DiveTypeName ) )```# Maps```{r}# check if background_ggoceanmaps existif (file.exists("../export/background_ggoceanmap.rds")) { trip <-readRDS("../export/background_ggoceanmap.rds")} else {# using ggOceanMaps trip <-basemap(limits =c(170, -110, 20, 59),bathymetry =TRUE,shapefiles ="Arctic",rotate =TRUE,grid.col =NA )# Make the graticules: lims <-attributes(trip)$limits graticule <- sf::st_graticule(c(lims[1], lims[3], lims[2], lims[4]),crs =attributes(trip)$proj,lon =seq(-180, 180, 45),lat =seq(-90, 90, 10) )# Plot trip = trip +geom_sf(data = graticule, color ="grey50") trip$layers <- trip$layers[c(1,3,2)]# save resultsaveRDS(trip, "../export/background_ggoceanmap.rds")}```::: panel-tabset## All```{r}# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat) & DiveTypeName =="Benthic") %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName)``````{r fig-map-all}#| fig-cap: "Kernel density plots of the all benthic dives performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park."#| fig-height: 5#| fig-asp: 0.6# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(x = df_kernel_dens_sf,# boundary = ne_coastline(scale = 110, returnclass = "sf"),H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plottest = trip +# new scalenew_scale_fill() +# kernel# ggplot(data=ne_coastline(scale = 110, returnclass = "sf"))+geom_sf()+geom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# kde_test,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# titlelabs(title ="All Benthic Dives",subtitle =paste(nrow(df_kernel_dens), "dives"))# reorder layerstest$layers = test$layers[c(1,2,4,3)]# plot test```## diff\<50```{r}# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat) & DiveTypeName_50 =="Benthic") %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName)``````{r fig-map-50}#| fig-cap: "Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 50 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park."#| fig-height: 5#| fig-asp: 0.6# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plottest = trip +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# titlelabs(title ="Benthic Dives with bathy-max_depth<50",subtitle =paste(nrow(df_kernel_dens), "dives"))# reorder layerstest$layers = test$layers[c(1,2,4,3)]# plot test```## diff\<100```{r}# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat) & DiveTypeName_100 =="Benthic") %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName)``````{r fig-map-100}#| fig-cap: "Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 100 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park."#| fig-height: 5#| fig-asp: 0.6# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plottest = trip +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# titlelabs(title ="Benthic Dives with bathy-max_depth<100",subtitle =paste(nrow(df_kernel_dens), "dives"))# reorder layerstest$layers = test$layers[c(1,2,4,3)]# plot test```## diff\<150```{r}# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat) & DiveTypeName_150 =="Benthic") %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName)``````{r fig-map-150}#| fig-cap: "Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 150 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park."#| fig-height: 5#| fig-asp: 0.6# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plottest = trip +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# titlelabs(title ="Benthic Dives with bathy-max_depth<150",subtitle =paste(nrow(df_kernel_dens), "dives"))# reorder layerstest$layers = test$layers[c(1,2,4,3)]# plot test```## diff\<200```{r}# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat) & DiveTypeName_200 =="Benthic") %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName)``````{r fig-map-200}#| fig-cap: "Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 200 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park."#| fig-height: 5#| fig-asp: 0.6# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plottest = trip +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# titlelabs(title ="Benthic Dives with bathy-max_depth<200",subtitle =paste(nrow(df_kernel_dens), "dives"))# reorder layerstest$layers = test$layers[c(1,2,4,3)]# plot test```:::# Windrose::: panel-tabset## All```{r fig-windrose-all}#| fig-cap: "Circular histogram plots displaying the time when they perform benthic dives across their foraging trip."#| # let's add a column with the local_timedata_windrose <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName =="Benthic")# figurewindrose_benthic <- data_windrose %>%group_by(id) %>%filter(DiveTypeName =="Benthic") %>%mutate(time =as.numeric(format(date_tz, format ="%H"))) %>%group_by(time) %>%summarise(nb_ind_hour =n()) %>%mutate(origin ="Benthic") %>%windrose(., grid =c(3000, 6000, 9000), legend_position ="top", facet = T)# displaywindrose_benthic``````{r}# perform testtest_benthic =r.test(conversion.circular(circular( data_windrose %>%group_by(id) %>%filter(DiveTypeName =="Benthic") %>%mutate(time =hour(date_tz) +minute(date_tz)/60+second(date_tz)/(60*60)) %>%pull(time),units ="hours" ),units ="radians"))# convert back into human reading hourconv_test =conversion.circular(circular(test_benthic$r.bar, units ="rad"), units ="hours")# printprint(paste0("The average time is ", round(conv_test[[1]],1), " (p-value: ", round(test_benthic$p.value, 1), ")"))```## Diff\<50```{r fig-windrose-50}#| fig-cap: "Circular histogram plots displaying the time when they perform benthic dives (diff<50) across their foraging trip."#| # let's add a column with the local_timedata_windrose <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_50 =="Benthic")# figurewindrose_benthic <- data_windrose %>%group_by(id) %>%filter(DiveTypeName =="Benthic") %>%mutate(time =as.numeric(format(date_tz, format ="%H"))) %>%group_by(time) %>%summarise(nb_ind_hour =n()) %>%mutate(origin ="Benthic") %>%windrose(., grid =c(1000, 2000, 3000), legend_position ="top", facet = T)# displaywindrose_benthic``````{r}# perform testtest_benthic =r.test(conversion.circular(circular( data_windrose %>%group_by(id) %>%filter(DiveTypeName_50 =="Benthic") %>%mutate(time =hour(date_tz) +minute(date_tz)/60+second(date_tz)/(60*60)) %>%pull(time),units ="hours" ),units ="radians"))# convert back into human reading hourconv_test =conversion.circular(circular(test_benthic$r.bar, units ="rad"), units ="hours")# printprint(paste0("The average time is ", round(conv_test[[1]],1), " (p-value: ", round(test_benthic$p.value, 1), ")"))```## Diff\<100```{r fig-windrose-100}#| fig-cap: "Circular histogram plots displaying the time when they perform benthic dives (diff<100) across their foraging trip."#| # let's add a column with the local_timedata_windrose <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_100 =="Benthic")# figurewindrose_benthic <- data_windrose %>%group_by(id) %>%filter(DiveTypeName =="Benthic") %>%mutate(time =as.numeric(format(date_tz, format ="%H"))) %>%group_by(time) %>%summarise(nb_ind_hour =n()) %>%mutate(origin ="Benthic") %>%windrose(., grid =c(1000, 2000, 3000), legend_position ="top", facet = T)# displaywindrose_benthic``````{r}# perform testtest_benthic =r.test(conversion.circular(circular( data_windrose %>%group_by(id) %>%filter(DiveTypeName_100 =="Benthic") %>%mutate(time =hour(date_tz) +minute(date_tz)/60+second(date_tz)/(60*60)) %>%pull(time),units ="hours" ),units ="radians"))# convert back into human reading hourconv_test =conversion.circular(circular(test_benthic$r.bar, units ="rad"), units ="hours")# printprint(paste0("The average time is ", round(conv_test[[1]],1), " (p-value: ", round(test_benthic$p.value, 1), ")"))```## Diff\<150```{r fig-windrose-150}#| fig-cap: "Circular histogram plots displaying the time when they perform benthic dives (diff<150) across their foraging trip."#| # let's add a column with the local_timedata_windrose <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_150 =="Benthic")# figurewindrose_benthic <- data_windrose %>%group_by(id) %>%filter(DiveTypeName =="Benthic") %>%mutate(time =as.numeric(format(date_tz, format ="%H"))) %>%group_by(time) %>%summarise(nb_ind_hour =n()) %>%mutate(origin ="Benthic") %>%windrose(., grid =c(1000, 2000, 3000), legend_position ="top", facet = T)# displaywindrose_benthic``````{r}# perform testtest_benthic =r.test(conversion.circular(circular( data_windrose %>%group_by(id) %>%filter(DiveTypeName_150 =="Benthic") %>%mutate(time =hour(date_tz) +minute(date_tz)/60+second(date_tz)/(60*60)) %>%pull(time),units ="hours" ),units ="radians"))# convert back into human reading hourconv_test =conversion.circular(circular(test_benthic$r.bar, units ="rad"), units ="hours")# printprint(paste0("The average time is ", round(conv_test[[1]],1), " (p-value: ", round(test_benthic$p.value, 1), ")"))```## Diff\<200```{r fig-windrose-200}#| fig-cap: "Circular histogram plots displaying the time when they perform benthic dives (diff<200) across their foraging trip."#| # let's add a column with the local_timedata_windrose <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_200 =="Benthic")# figurewindrose_benthic <- data_windrose %>%group_by(id) %>%filter(DiveTypeName =="Benthic") %>%mutate(time =as.numeric(format(date_tz, format ="%H"))) %>%group_by(time) %>%summarise(nb_ind_hour =n()) %>%mutate(origin ="Benthic") %>%windrose(., grid =c(1000, 2000, 3000), legend_position ="top", facet = T)# displaywindrose_benthic``````{r}# perform testtest_benthic =r.test(conversion.circular(circular( data_windrose %>%group_by(id) %>%filter(DiveTypeName_200 =="Benthic") %>%mutate(time =hour(date_tz) +minute(date_tz)/60+second(date_tz)/(60*60)) %>%pull(time),units ="hours" ),units ="radians"))# convert back into human reading hourconv_test =conversion.circular(circular(test_benthic$r.bar, units ="rad"), units ="hours")# printprint(paste0("The average time is ", round(conv_test[[1]],1), " (p-value: ", round(test_benthic$p.value, 1), ")"))```:::# Boxplot of the benthic dive rate (day vs. night)Since we thought of accoutering for the variation in dive duration to explore how benthic dive rate varies across day-time vs night-time, I just looked at the difference in dive duration between both period of the day.```{r}data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%ggplot(aes(x=day_night, y=Dduration)) +# the boxplotgeom_boxplot(outlier.shape =NA) +# zoomcoord_cartesian(ylim =c(0, 2600)) +# labelslabs(x ="Time of the day", y ="Dive duration (s)")``````{r}broom::tidy(wilcox.test( data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# at nightfilter(day_night =="night") %>%# get dive durationpull(Dduration), data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# at nightfilter(day_night =="day") %>%# get dive durationpull(Dduration)))```::: panel-tabset## All### Graph {-}```{r fig-boxplot-all}#| fig-cap: "Boxplot of the number of (all) benthic dives per hour and individual accross day-time and night-time"#| data_plot = data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName =="Benthic") %>%mutate(hour = lubridate::hour(date),day =as.Date(date)) %>%group_by(id, day, hour, day_night) %>%summarise(nb_benthic_dive =n(), mean_dive_duration =mean(Dduration),.groups ="drop") %>%group_by(id, day_night) %>%summarise(nb_average_benthic_dive =mean(nb_benthic_dive),mean_dive_duration =mean(mean_dive_duration),.groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +geom_boxplot(outlier.shape =NA) +labs(x ="Time of the day", y ="Average number of benthic dives per seal per hour")```### Model {-}```{r fig-fit-data-all}#| fig-cap: Result from `descdist` function to determine the best distribution to use for data modelling# determine what distribution would best describe our datadescdist(data_plot$nb_average_benthic_dive, discrete =FALSE)```Based on @fig-fit-data-all, it seems that a `gamma` distribution would be the best to model this dataset.$$ HourlyBenthicDiveRate \sim DayNight + MeanDiveDuration + (1|ID)$$ {#eq-mdl-benthic-rate}Based on a table shape like the following:| ID | HourlyBenthicDiveRate | DayNight | MeanDiveDuration ||-----|-----------------|----------|------------------|| A | 10 | Night | 20 || A | 5 | Day | 21 || B | 11 | Night | 20 || B | 4 | Day | 22 |```{r}# modelmdl_all =glmer(nb_average_benthic_dive ~ day_night +scale(mean_dive_duration) + (1|id), family = Gamma,data = data_plot %>%mutate(day_night =factor(day_night,# ordered = T,levels=c("night","day"))))# summarybroom.mixed::tidy(mdl_all) %>%gt() %>%tab_header(title ="GLMM")```### Go Further... {-}```{r}data_plot = data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName =="Benthic") %>%# I'm not so sure about this part..mutate(period_duration =if_else( day_night =="night",getNightLength(date, Lat, Long),24-getNightLength(date, Lat, Long) )) %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(),mean_dive_duration =mean(Dduration),mean_period_duration =mean(period_duration),.groups ="drop" )```$$ NbBenthicDive \sim DayNight + MeanDiveDuration + MeanDayNightDuration + (1|ID)$$ Based on a table shape like the following:| ID | NbBenthic | DayNight | MeanDiveDuration | MeanDayNightDuration||-----|-----------------|----------|------------------|---------------|| A | 10 | Night | 20 |7 || A | 5 | Day | 21 |17 || B | 11 | Night | 20 |6 || B | 4 | Day | 22 |18 |For more information: <https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#overdispersion>```{r}# modelmdl_all_extra =glmer(nb_benthic_dive ~ day_night +scale(mean_dive_duration) +scale(mean_period_duration) + (1|id), family ="poisson",data = data_plot %>%mutate(day_night =factor(day_night,# ordered = T,levels=c("day","night"))))# summarybroom.mixed::tidy(mdl_all_extra) %>%gt() %>%tab_header(title ="GLMM")``````{r}qqnorm(residuals(mdl_all_extra))qqline(residuals(mdl_all_extra))```## diff\<50### Graph {-}```{r fig-boxplot-50}#| fig-cap: "Boxplot of the number of (diff<50) benthic dives per hour and individual accross day-time and night-time"#| data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_50 =="Benthic") %>%mutate(hour = lubridate::hour(date),day =as.Date(date)) %>%group_by(id, day, hour, day_night) %>%summarise(nb_benthic_dive =n(), mean_dive_duration =mean(Dduration),.groups ="drop") %>%group_by(id, day_night) %>%summarise(nb_average_benthic_dive =mean(nb_benthic_dive),mean_dive_duration =mean(mean_dive_duration),.groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +geom_boxplot(outlier.shape =NA) +labs(x ="Time of the day", y ="Average number of benthic dives per seal per hour")```### Model {-}```{r fig-fit-data-50}#| fig-cap: Result from `descdist` function to determine the best distribution to use for data modelling# determine what distribution would best describe our datadescdist(data_plot$nb_average_benthic_dive, discrete =FALSE)```Based on @fig-fit-data-50, it seems that a `gamma` distribution would be the best to model this dataset.$$ HourlyBenthicDiveRate \sim DayNight + MeanDiveDuration + (1|ID)$$ {#eq-mdl-benthic-rate}Based on a table shape like the following:| ID | HourlyBenthicDiveRate | DayNight | MeanDiveDuration ||-----|-----------------|----------|------------------|| A | 10 | Night | 20 || A | 5 | Day | 21 || B | 11 | Night | 20 || B | 4 | Day | 22 |```{r}# modelmdl_50 =glmer(nb_average_benthic_dive ~ day_night +scale(mean_dive_duration) + (1|id), family = Gamma,data = data_plot %>%mutate(day_night =factor(day_night,levels=c("night","day"))))# summarybroom.mixed::tidy(mdl_50) %>%gt() %>%tab_header(title ="GLMM")```### Go Further... {-}```{r}data_plot = data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_50 =="Benthic") %>%# I'm not so sure about this part..mutate(period_duration =if_else( day_night =="night",getNightLength(date, Lat, Long),24-getNightLength(date, Lat, Long) )) %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(),mean_dive_duration =mean(Dduration),mean_period_duration =mean(period_duration),.groups ="drop" )```For more information: <https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#overdispersion>$$ NbBenthicDive \sim DayNight + MeanDiveDuration + MeanDayNightDuration + (1|ID)$$ Based on a table shape like the following:| ID | NbBenthic | DayNight | MeanDiveDuration | MeanDayNightDuration||-----|-----------------|----------|------------------|---------------|| A | 10 | Night | 20 |7 || A | 5 | Day | 21 |17 || B | 11 | Night | 20 |6 || B | 4 | Day | 22 |18 |```{r}# modelmdl_50_extra =glmer(nb_benthic_dive ~ day_night +scale(mean_dive_duration) +scale(mean_period_duration) + (1|id), family ="poisson",data = data_plot %>%mutate(day_night =factor(day_night,# ordered = T,levels=c("day","night"))))# summarybroom.mixed::tidy(mdl_50_extra) %>%gt() %>%tab_header(title ="GLMM")``````{r}qqnorm(residuals(mdl_50_extra))qqline(residuals(mdl_50_extra))```## diff\<100### Graph {-}```{r fig-boxplot-100}#| fig-cap: "Boxplot of the number of (diff<100) benthic dives per hour and individual accross day-time and night-time"#| data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_100 =="Benthic") %>%mutate(hour = lubridate::hour(date),day =as.Date(date)) %>%group_by(id, day, hour, day_night) %>%summarise(nb_benthic_dive =n(), mean_dive_duration =mean(Dduration),.groups ="drop") %>%group_by(id, day_night) %>%summarise(nb_average_benthic_dive =mean(nb_benthic_dive),mean_dive_duration =mean(mean_dive_duration),.groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +geom_boxplot(outlier.shape =NA) +labs(x ="Time of the day", y ="Average number of benthic dives per seal per hour")```### Model {-}```{r fig-fit-data-100}#| fig-cap: Result from `descdist` function to determine the best distribution to use for data modelling# determine what distribution would best describe our datadescdist(data_plot$nb_average_benthic_dive, discrete =FALSE)```Based on @fig-fit-data-100, it seems that a `gamma` distribution would be the best to model this dataset.$$ HourlyBenthicDiveRate \sim DayNight + MeanDiveDuration + (1|ID)$$ {#eq-mdl-benthic-rate}Based on a table shape like the following:| ID | HourlyBenthicDiveRate | DayNight | MeanDiveDuration ||-----|-----------------|----------|------------------|| A | 10 | Night | 20 || A | 5 | Day | 21 || B | 11 | Night | 20 || B | 4 | Day | 22 |```{r}# modelmdl_100 =glmer(nb_average_benthic_dive ~ day_night +scale(mean_dive_duration) + (1|id), family = Gamma,data = data_plot %>%mutate(day_night =factor(day_night,levels=c("night","day"))))# summarybroom.mixed::tidy(mdl_100) %>%gt() %>%tab_header(title ="GLMM")```### Go Further... {-}```{r}data_plot = data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_100 =="Benthic") %>%# I'm not so sure about this part..mutate(period_duration =if_else( day_night =="night",getNightLength(date, Lat, Long),24-getNightLength(date, Lat, Long) )) %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(),mean_dive_duration =mean(Dduration),mean_period_duration =mean(period_duration),.groups ="drop" )```For more information: <https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#overdispersion>$$ NbBenthicDive \sim DayNight + MeanDiveDuration + MeanDayNightDuration + (1|ID)$$ Based on a table shape like the following:| ID | NbBenthic | DayNight | MeanDiveDuration | MeanDayNightDuration||-----|-----------------|----------|------------------|---------------|| A | 10 | Night | 20 |7 || A | 5 | Day | 21 |17 || B | 11 | Night | 20 |6 || B | 4 | Day | 22 |18 |```{r}# modelmdl_100_extra =glmer(nb_benthic_dive ~ day_night +scale(mean_dive_duration) +scale(mean_period_duration) + (1|id), family ="poisson",data = data_plot %>%mutate(day_night =factor(day_night,# ordered = T,levels=c("day","night"))))# summarybroom.mixed::tidy(mdl_50_extra) %>%gt() %>%tab_header(title ="GLMM")``````{r}qqnorm(residuals(mdl_100_extra))qqline(residuals(mdl_100_extra))```## diff\<150### Graph {-}```{r fig-boxplot-150}#| fig-cap: "Boxplot of the number of (diff<150) benthic dives per hour and individual accross day-time and night-time"#| data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_150 =="Benthic") %>%mutate(hour = lubridate::hour(date),day =as.Date(date) ) %>%group_by(id, day, hour, day_night) %>%summarise(nb_benthic_dive =n(),mean_dive_duration =mean(Dduration),.groups ="drop" ) %>%group_by(id, day_night) %>%summarise(nb_average_benthic_dive =mean(nb_benthic_dive),mean_dive_duration =mean(mean_dive_duration),.groups ="drop" )data_plot %>%ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +geom_boxplot(outlier.shape =NA) +labs(x ="Time of the day",y ="Average number of benthic dives per seal per hour" )```### Model {-}```{r fig-fit-data-150}#| fig-cap: Result from `descdist` function to determine the best distribution to use for data modelling# determine what distribution would best describe our datadescdist(data_plot$nb_average_benthic_dive, discrete =FALSE)```Based on @fig-fit-data-150, it seems that a `gamma` distribution would be the best to model this dataset.$$ HourlyBenthicDiveRate \sim DayNight + MeanDiveDuration + (1|ID)$$ {#eq-mdl-benthic-rate}Based on a table shape like the following:| ID | HourlyBenthicDiveRate | DayNight | MeanDiveDuration ||-----|-----------------|----------|------------------|| A | 10 | Night | 20 || A | 5 | Day | 21 || B | 11 | Night | 20 || B | 4 | Day | 22 |```{r}# modelmdl_150 =glmer(nb_average_benthic_dive ~ day_night +scale(mean_dive_duration) + (1|id), family = Gamma,data = data_plot %>%mutate(day_night =factor(day_night,levels=c("night","day"))))# summarybroom.mixed::tidy(mdl_150) %>%gt() %>%tab_header(title ="GLMM")```### Go Further... {-}```{r}data_plot = data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_150 =="Benthic") %>%# I'm not so sure about this part..mutate(period_duration =if_else( day_night =="night",getNightLength(date, Lat, Long),24-getNightLength(date, Lat, Long) )) %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(),mean_dive_duration =mean(Dduration),mean_period_duration =mean(period_duration),.groups ="drop" )```For more information: <https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#overdispersion>$$ NbBenthicDive \sim DayNight + MeanDiveDuration + MeanDayNightDuration + (1|ID)$$ Based on a table shape like the following:| ID | NbBenthic | DayNight | MeanDiveDuration | MeanDayNightDuration||-----|-----------------|----------|------------------|---------------|| A | 10 | Night | 20 |7 || A | 5 | Day | 21 |17 || B | 11 | Night | 20 |6 || B | 4 | Day | 22 |18 |```{r}# modelmdl_150_extra =glmer(nb_benthic_dive ~ day_night +scale(mean_dive_duration) +scale(mean_period_duration) + (1|id), family ="poisson",data = data_plot %>%mutate(day_night =factor(day_night,# ordered = T,levels=c("day","night"))))# summarybroom.mixed::tidy(mdl_150_extra) %>%gt() %>%tab_header(title ="GLMM")``````{r}qqnorm(residuals(mdl_150_extra))qqline(residuals(mdl_150_extra))```## diff\<200### Graph {-}```{r fig-boxplot-200}#| fig-cap: "Boxplot of the number of (diff<200) benthic dives per hour and individual accross day-time and night-time"#| data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_200 =="Benthic") %>%mutate(hour = lubridate::hour(date),day =as.Date(date)) %>%group_by(id, day, hour, day_night) %>%summarise(nb_benthic_dive =n(), mean_dive_duration =mean(Dduration),.groups ="drop") %>%group_by(id, day_night) %>%summarise(nb_average_benthic_dive =mean(nb_benthic_dive),mean_dive_duration =mean(mean_dive_duration),.groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +geom_boxplot(outlier.shape =NA) +labs(x ="Time of the day", y ="Average number of benthic dives per seal per hour")```### Model {-}```{r fig-fit-data-200}#| fig-cap: Result from `descdist` function to determine the best distribution to use for data modelling# determine what distribution would best describe our datadescdist(data_plot$nb_average_benthic_dive, discrete =FALSE)```Based on @fig-fit-data-200, it seems that a `gamma` distribution would be the best to model this dataset.$$ HourlyBenthicDiveRate \sim DayNight + MeanDiveDuration + (1|ID)$$ {#eq-mdl-benthic-rate}Based on a table shape like the following:| ID | HourlyBenthicDiveRate | DayNight | MeanDiveDuration ||-----|-----------------|----------|------------------|| A | 10 | Night | 20 || A | 5 | Day | 21 || B | 11 | Night | 20 || B | 4 | Day | 22 |```{r}# modelmdl_200 =glmer(nb_average_benthic_dive ~ day_night +scale(mean_dive_duration) + (1|id), family = Gamma,data = data_plot %>%mutate(day_night =factor(day_night,levels=c("night","day"))))# summarybroom.mixed::tidy(mdl_200) %>%gt() %>%tab_header(title ="GLMM")```### Go Further... {-}```{r}data_plot = data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_200 =="Benthic") %>%# I'm not so sure about this part..mutate(period_duration =if_else( day_night =="night",getNightLength(date, Lat, Long),24-getNightLength(date, Lat, Long) )) %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(),mean_dive_duration =mean(Dduration),mean_period_duration =mean(period_duration),.groups ="drop" )```For more information: <https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#overdispersion>$$ NbBenthicDive \sim DayNight + MeanDiveDuration + MeanDayNightDuration + (1|ID)$$ Based on a table shape like the following:| ID | NbBenthic | DayNight | MeanDiveDuration | MeanDayNightDuration||-----|-----------------|----------|------------------|---------------|| A | 10 | Night | 20 |7 || A | 5 | Day | 21 |17 || B | 11 | Night | 20 |6 || B | 4 | Day | 22 |18 |```{r}# modelmdl_200_extra =glmer(nb_benthic_dive ~ day_night +scale(mean_dive_duration) +scale(mean_period_duration) + (1|id), family ="poisson",data = data_plot %>%mutate(day_night =factor(day_night,# ordered = T,levels=c("day","night"))))# summarybroom.mixed::tidy(mdl_200_extra) %>%gt() %>%tab_header(title ="GLMM")``````{r}qqnorm(residuals(mdl_200_extra))qqline(residuals(mdl_200_extra))```:::# Boxplot of the number of benthic dive (day vs. night)::: panel-tabset## All```{r fig-boxplotbis-all}#| fig-cap: "Boxplot of the average number of (all) benthic dives per hour and individual accross day-time and night-time"#| data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName =="Benthic") %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_benthic_dive)) +geom_boxplot(outlier.shape =NA) +coord_cartesian(ylim=c(0,1000))+labs(x ="Time of the day", y ="Number of benthic dives per individual")``````{r}broom::tidy(glm( nb_benthic_dive ~ day_night,family ="quasipoisson",data = data_plot)) %>%gt() %>%fmt_number(columns =2:5,decimals =2,suffixing =TRUE) %>%tab_header(title ="nb_benthic_dive ~ day_night")```## diff\<50```{r fig-boxplotbis-50}#| fig-cap: "Boxplot of the average number of (diff<50) benthic dives per hour and individual accross day-time and night-time"#| data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_50 =="Benthic") %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_benthic_dive)) +geom_boxplot(outlier.shape =NA) +coord_cartesian(ylim=c(0,250))+labs(x ="Time of the day", y ="Number of benthic dives per individual")``````{r}broom::tidy(glm( nb_benthic_dive ~ day_night,family ="quasipoisson",data = data_plot)) %>%gt() %>%fmt_number(columns =2:5,decimals =2,suffixing =TRUE) %>%tab_header(title ="nb_benthic_dive ~ day_night")```## diff\<100```{r fig-boxplotbis-100}#| fig-cap: "Boxplot of the average number of (diff<100) benthic dives per hour and individual accross day-time and night-time"#| data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_100 =="Benthic") %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_benthic_dive)) +geom_boxplot(outlier.shape =NA) +coord_cartesian(ylim=c(0,250))+labs(x ="Time of the day", y ="Number of benthic dives per individual")``````{r}broom::tidy(glm( nb_benthic_dive ~ day_night,family ="quasipoisson",data = data_plot)) %>%gt() %>%fmt_number(columns =2:5,decimals =2,suffixing =TRUE) %>%tab_header(title ="nb_benthic_dive ~ day_night")```## diff\<150```{r fig-boxplotbis-150}#| fig-cap: "Boxplot of the average number of (diff<150) benthic dives per hour and individual accross day-time and night-time"#| data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_150 =="Benthic") %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_benthic_dive)) +geom_boxplot(outlier.shape =NA) +coord_cartesian(ylim=c(0,250))+labs(x ="Time of the day", y ="Number of benthic dives per individual")``````{r}broom::tidy(glm( nb_benthic_dive ~ day_night,family ="quasipoisson",data = data_plot)) %>%gt() %>%fmt_number(columns =2:5,decimals =2,suffixing =TRUE) %>%tab_header(title ="nb_benthic_dive ~ day_night")```## diff\<200```{r fig-boxplotbis-200}#| fig-cap: "Boxplot of the average number of (diff<200) benthic dives per hour and individual accross day-time and night-time"#| data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_200 =="Benthic") %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_benthic_dive)) +geom_boxplot(outlier.shape =NA) +coord_cartesian(ylim=c(0,250))+labs(x ="Time of the day", y ="Number of benthic dives per individual")``````{r}broom::tidy(glm( nb_benthic_dive ~ day_night,family ="quasipoisson",data = data_plot)) %>%gt() %>%fmt_number(columns =2:5,decimals =2,suffixing =TRUE) %>%tab_header(title ="nb_benthic_dive ~ day_night")```:::# ExtraLet's look at the distribution of the dive duration to decide a threshold.```{r}data_dive %>%ggplot(aes(x = Dduration)) +geom_histogram(fill ="grey80",color ="black") +geom_vline(xintercept =23*60, linetype="dashed", col ="lightgreen") +geom_vline(xintercept =4000, linetype="dashed", col ="salmon") +labs(y ="# of dives", x ="Dive duration (s)") +theme_minimal() +theme(axis.line =element_line(arrow =arrow(type ="closed", length =unit(0.1, "inches"))))```Ok, so let's pick 2500 sec and look on a map where this dives are.```{r}# define palettedata_inter = data_dive %>%filter(!is.na(Lat) & Dduration >4000) %>%arrange(id, date) %>%mutate(date_post = date + Dduration,diff_time =as.numeric(c(difftime(tail(date, -1), head(date_post, -1)),0), units ="mins"))colPal <- scales::col_numeric(palette ="RdYlGn",domain =c(0, max(data_inter$diff_time)))trip +# new scalenew_scale_fill() +# kernel ggspatial::geom_spatial_point(data = data_inter,aes(x = Long, y = Lat),crs =4326,col =colPal(data_inter$diff_time),size =1 ) +# no display alphaguides() +# titlelabs(title ="Dives longer than 4000 seconds",subtitle =paste0(nrow(data_inter), " dives"))```